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Abstract 

°t7Results  of  (i)  a  numerical  investigation,  based  on  an  energy  consistent 
moving-singularity  dynamic  finite  element  procedure,  of  fast  crack  propagation 
in  a  finite  plate;  (ii)  numerical  simulation  of  experimental  data  on  fast 
crack  propagation  and  arrest  in  a  double-cantilever-beam  specimen;  and  (iii) 
stress-intensity  factor  solutions  in  a  thermally  shocked  cylindrical  vessel 


containing  an  inner  surface  (meridional)  elliptical  flaw,  are  presented. 
Comparison  of  these  results  with  other  available  solutions,  and  pertinent 
discussions,  are  included. 

Introduction 

Concise  summaries  of  the  current  status  of  the  subject  of  dynamic  crack 

propagation  can  be  found  in  recent  review  articles  by  Achenbach  [1]  and  Freund 

[2],  In  Refs.  1  and  2  several  analytical  solutions  of  linear  elasto-dynamic 

equations  for  crack  propagation  in  plane  bodies  with  infinite  domains  have 

been  reviewed.  For  finite  bodies  containing  cracks  and  subjected  to  time- 

dependent  loading,  the  interactior  with  a  crack-tip  of  stress-waves  reflected 

from  the  boundaries  and/or  emanated  by  the  other  moving  crack-tip  play  an 

important  role  in  determing  the  intensity  of  the  dynamic  singular  stress-field 

at  the  considered  crack-tip.  Because  of  the  analytical  intractability  of 

such  elasto-dynamic  problems  for  finite  domains,  computational  techniques  are 
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mandatory.  A  critical  appraisal  of  several  different  numerical  (finite  element 
and  finite  difference)  techniques  was  made  by  Kanninen  [3]  in  1978.  Most  of 
the  finite  element  methods  reviewed  in  [3]  use  conventional  assumed-displacement 
finite  elements  near  the  crack-tip  and  hence  do  not  accout  for  the  known 
crack-tip  singularity.  Moreover  in  these  methods,  crack-propagation  is 
simulated  by  the  well-known  "node-release"  technique,  which  was  argued  in 
[3]  to  be  not  sufficiently  accurate.  The  literature  an  dynamic  finite  element 
methods  for  simulation  of  fast  fracture,  since  the  appearance  of  [3],  has 
been  reviewed  in  [4,5,6]. 

In  Refs.  [4,5,6]  the  authors  have  presented  a  "moving-singularity" 
finite  element  procedure  for  the  analysis  of  dynamic  crack  propagation  in 
arbitrarily  shaped  finite  bodies.  In  this  procedure  a  singular  crack-tip 
element,  within  which  a  large  number  of  analytical  eigen-functions  correspond¬ 
ing  to  a  propagating  crack  [4,7]  are  used  as  basis  functions  for  displacements, 
translates  by  an  arbitrary  amount  A T.  in  each  time-increment  At  of  the  numerical 
time-integration  procedure.  The  moving  crack-tip  singular-element,  within 
which  the  crack-tip  always  has  a  fixed  location,  retains  its  shape  at  all 
times,  but  the  mesh  of  conventional  (isoparametric)  finite  elements,  surrounding 
the  moving  singular- element ,  deforms  accordingly.  An  energy-consistent 
variational  statement  was  developed  in  [4,5,6]  as  a  basis  for  the  above  moving 
singularity  procedure.  It  was  also  shown  in  [4,5,6]  that  the  procedure  there 
in  lead  to  a  direct  evaluation  of  the  dynamic  stress-intensity  factors  for 
propagating  cracks.  Several  numerical  studies  were  presented  in  [4,6]  to 
illustrate  the  relative  efficiency  of  the  above  procedure  as  compared  to  the 
"node-release"  techniques  reviewed  in  [3,4]. 

In  the  present  paper  we  present  further  numerical  results  for  the 
problem  of  dynamic  propagation,  at  different  constant  velocities,  of  a  cen¬ 
trally  located  crack  in  a  square  panel.  These  results  augment  the  pre- 


liminary  conclusions  reached  in  [4,6]  concerning  the  effects  of  interactions 

of  stress-waves  emanated  from  a  moving  crack-tip  and  those  reflected  from  the 
boundaries  of  the  panel  on  the  dynamic  stress-intensity  factor  at  the  con¬ 
sidered  crack-tip.  These  results  are  also  shown  to  lead  to  a  simple  formula 
for  estimating  the  dynamic  K  factor  in  similar  situations.  Also  presented 
are  the  results  of  simulation  of  data  on  the  crack-tip  velocity  history  in  an 
experiment  on  a  double-cantilever-beam  specimen  reported  by  Kalthoff  et  al 
[8].  The  dynamic  stress-intensity  factors,  for  this  specimen,  computed  from 
the  present  procedure  are  compared  with  the  experimentally  determined  data  of 
[8],  and  independent  numerical  results  of  Kobayashi  et  al  [9]  and  Popelar  et 
al  [10],  and  pertinent  discussions  are  presented.  j- 

In  the  second  part  of  the  paper,  results  of  an  investigation  of  the 
stress-intensity  factors  near  the  border  of  a  meridional  semi-elliptical 
surface  flaw  at  the  inner  surface  of  a  cylindrical  pressure  vessel  which  is 
subjected  to  a  thermal  shock,  are  presented.  The  presently  reported  results 
are  obtained  by  using  "three-dimensional  hybrid  crack-elements"  near  the 
crack  front,  the  development  of  which  was  reported  eariler  by  the  authors 
[11,12,13].  The  present  results  are  compared  with  those  reported  earlier 
by  Kobayashi  et  al  [14]. 


Part  I.  Dynamic  Crack  Propagation  Analysis 

Synopsis  of  the  Analysis  Procedure:  In  the  procedure  adopted  in  the 
paper,  the  basis  functions  used  for  displacement,  velocity,  and  acceleration 
in  the  crack-tip  "singular  element"  are: 


u  (r,x  ,t)  =  u  ,(f,x„,v)g  (t)  la=l,2  ;  j=l,...N] 
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where  u^.  correspond  to  the  "steady-state"  (ie.,  which  are  invariant  to  an 
observer  moving  with  the  crack-tip)  eigen-function  solutions  for  the  elasto-- 
dynamic  wave  equations  (with  independent  variables  £,  and  X2)  for  crack  propaga¬ 
tion  at  constant  velocity  v  in  a  plane  domain.  Note  that  xa(a=l,2)  are  fixed 

coordinates,  with  X2  =  0  defining  the  crack  plane  and  £  =  x^  -  vt.  It  is  noted 

-1/2 

that  the  first  term,  viz. ,  u  ,  leads  to  the  appropriate  (r  )  type  singularity 
in  strains  and  stresses.  The  singular  element  in  the  present  procedure  is 
surrounded  by  the  usual  isoparametric  [8-noded,  in  the  present  case]  elements. 

The  displacement  compatibility  between  the  singular  element  and  the  surrounding 
isoparametric  elements  is  satisfied  in  the  present  analysis  through  a  least- 
square  technique. 

Consider  two  instants  of  time,  t^  and  t2  *  +  At.  Assume  that  in  a 

mode  I  crack  propagation  problem,  the  crack- lengths  ar  t^  and  t2  are,  respec¬ 
tively,  E^  and  E^  +  AE.  Let  the  displacements,  strains,  and  stresses  at  t^ 
be  denoted  by  ul,  >  and  °ij ’  resPectively »  while  those  at  t2  are  denoted 
by  a  superscript  two  for  each  variable.  The  variables  at  time  t^  are  pressumed 
to  be  known.  It  has  been  shown  in  [4,5]  that  the  variational  principle 
governing  the  dynamic  crack  propagation  between  times  t^  and  t?  can  be  written  as: 


+  a]  J5e2  +  pCu}  +  u2)6u2  dV 
ij  ij  ill 
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where  M  is  the  domain  of  the  body,  and  S  is  the  boundary  of  V~  where  tractions 
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are  prescribed,  at  time  t2*  are  the  prescribed  tractions  at  time  t^  at 

_2 

S  (~S  )  and  T.  are  the  prescribed  tractions  at  S  as  well  as  at  the  newly 
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created  crack  surface  AZ  at  time  t2*  It  is  seen  that  at  AZ  are  the 

cohesive  forces  holding  the  crack-faces  together  at  time  t^.  In  the  above, 
mode  I  conditions  are  assumed;  hence,  only  the  upper  half  of  the  domain  with 
the  crack  face  Z+  is  considered. 

In  the  variational  principle  in  Eq.  (4),  the  varialbes  u^,  and  are 

i  ij 

2  2  2  2 

pressumed  known;  while  cf..,  e,.,  and  u.  are  the  variables.  The  variables  u 

ij  ij  3  j 

are  assumed  according  to  Eq.  (1),  with  the  velocity  v2  appearing  in  them. 

Further,  the  variational  principle  in  Eq.  (4)  is  used  to  develop  a  discrete 

(finite  element)  approximation  for  a  (finite  element)  mesh  at  time  t^.  Note 

that  at  time  t2»  in  the  present  problem,  the  crack-tip  is  located  at  x  =  Z^ 

+  AZ  and  hence  the  present  crack-element  is  centered  at  =  Z^  +  AZ.  In 

developing  the  equations  for  the  finite  element  mesh  at  t2>  it  is  seen  from 

Eq.  (4),  that  the  variation  of  a'!',  and  u^  must  be  known  in  the  finite  element 

13  3 

mesh  at  time  t„.  However,  a . .  and  u. ,  and  were  solved  for,  in  the  finite 
2  13  3  3 

element  mesh  at  t^.  In  the  mesh  at  t^  the  crack-tip  is  located  at  =  Z^  and 
hence  the  crack-element  is  centered  at  Z^.  Thus,  between  t^  and  t2  (=  t^  + 

At)  the  crack-eiement  is  translated  by  an  amount  AZ.  While  the  crack-element 
is  translated,  only  the  elements  immediately  surrounding  the  moving  crack-tip 
are  distorted.  Thus  the  finite  element  meshes  at  times  t^  and  t2  differ  only 
in  the  location  of  crack-tip  (and  hence  the  crack-element)  and  the  shapes 
of  the  immediately  surrounding  isoparametric  elements.  Thus,  the  known  data 
for  olj  and  uj  in  the  mesh  at  t^  is  interpolated  easily  into  corresponding 
data  in  the  mesh  at  t2>  Based  on  these  concepts,  the  development  of  the  finite- 
element  equations  from  the  principle  in  Eq .  (4),  and  the  numerical  integration 


of  these  equations  follows  the  well-established  procedures.  Further  details 


can  be  found  in  [4,5]  where  it  is  shown  that  the  dynamic  k-factors  can  be 
computed  directly  in  the  present  analysis  procedure. 

Propagation  of  a  Central  Crack  in  a  Square  Panel 

In  Ref.  [6],  the  problem  of  a  centrally  cracked  square  panel  (L=W=40mm) 

10  2 

[with  material  properties:  u(shear  modulus)  =  2.94  x  10  N/mm  ;  v(Poisson 

3  3 

ratio)  =  0.286;  and  p(mass  density)  =  2.45  x  10  Kg/ra  ],  which  was  subject  to 
time- independent  tensile  stress  at  the  edges  of  the  specimen  parrellel  to  the 
crack  axis,  was  considered.  The  crack  was  assumed  to  start  propagating  from 
an  initial  length,  ZQ/W  =0.2,  and  to  grow  symmetrically  with  a  constant 
velocity  v.  Four  different  cases  of  v,  namely,  (v/Cg)  =  0.2,  0.4,  0.6  and 
0.8,  respectively  (where  Cg  is  the  shear  wave  speed)  were  considered.  These 
limited  results  appeared  to  suggest  a  simple  formula  to  estimate  the  dynamic 
stress-intensity  factor  in  such  problems.  Until  the  time  (t=Rc)  taken  by 
the  Rayliegh  waves  emanating  one  crack-tip  to  interact  with  the  other  moving 
crack-tip,  the  dynamic  stress  intensity  factor  was  noted  in  [6]  to  be  given 
approximately  by  the  equation: 

*oo 

K,  =  F(Z  )  K  (Z)  k(v)  t  <  R  (5) 

do  c 

where  F(Zq)  is  the  finite  size  correction  factor  in  the  static  stress  intensity 

*oo 

factor  K  for  the  given  geometry  and  loading  at  a  crack  length  Z  ;  K  is 
s  o 

the  static  factor  (which,  in  general,  is  not  equal  to  the  static  stress-intensity 

OO 

factor  K  )  for  an  infinite  body  subjected  to  uniform  stress  normal  to  crack 
s 

*oo 

axis;  and  k(v)  is  a  universal  velocity  factor.  The  expressions  for  K  and 

k(v)  were  given  by  Eshelby  [15],  and  Broberg  [16]  respectively.  After  the  time 

t  =  R  ,  the  dynamic  stress-intensity  factor  was  found  in  [6]  to  be  given  by 
c 


the  approximate  relation: 


(6) 


Kd  =  S  ^  (£)  K(v)  t  >  Rc 

*0O  f 

where  0  =  [F(£  )  K  (£)]/[K  (£)]  at  £  =  £„r 

O  s 

where  is  the  static  stress-intensity  factor  for  the  present  finite  domain 
s 

and  £  is  the  current  crack  length  at  t  =  R  .  It  is  seen  that 

£or  =  £  +  v  (2£  ) / (CR-v)  (7) 

RC  O  O  K 

where  C  is  the  Rayleigh  wave  speed. 

K 

Here  we  present  further  results  for  the  case  when  a  central  crack  in  a 
square  panel  starts  from  an  initial  length  of  (E^/W)  =  0.1  and  propagates  with 
a  constant  velocity  v.  Two  cases  of  v,  namely,  (v/Cg)  =  0.1  and  0.2,  respective¬ 
ly,  are  considered. 

The  finite  element  mesh  at  the  initial  crack  length  £q,  for  a  quadrant 
of  the  panel  (which  only  is  modelled  due  to  symmetry),  is  shown  in  Fig.  1. 

This  mesh  has  262  nodes  and  564  degrees  of  freedom  before  the  imposition  of 
the  appropriate  symmetry  conditions  at  the  boundaries.  The  computed  results 
for  the  normalized  dynamic  stress-intensity  factor  (normalized  by  a/fr£,  £ 
being  the  current  crack  length)  are  shown  in  Fig.  2. 

For  the  present  material,  (CR/Cg)  =  0.9238.  Thus  it  is  seen  from  Eq. 

(7)  that  for  the  present  cases  of  (£q/W)  =  0.1,  and  (v/Cg)  =0.1  and  0.2,  one 
obtains  (£Rc/W)  .12,  and  .15,  respectively.  Shown  in  Fig.  (2)  are  the 

curves:  (i)  the  normalized  static-stress-intensity  factor  ;  (ii)  the 

presently  computed  normalized  dynamic  stress-intensity  factor;  and  (iii)  the 
approximations  for  the  dynamic-stress  intensity  factor  as  given  by  Eqs.  (*') 
and  (7). 

It  is  seen  that  the  correlation  between  the  directly  computed  dynamic 
K-factor  and  that  obtained  from  the  simple  approximations  in  Eqs.  (M  and  (7) 
is  very  good.  Approximate  formulae  similar  to  those  given  in  Eq.  (6)  and  (7) 


for  the  more  common  dynamic  fracture  test  specimens  such  as  the  double 


cantilever  specimen,  wedge-loaded  single  edge  notch  specimen,  and  modified 
compact  tension  specimen  would  greatly  aid  in  the  estimation  of  dynamic  K- 
factors  from  the  experimentally  obtained  data  for  E(t)  and/or  dE/dt.  Such 
studies  are  currently  underway  and  will  be  reported  elsewhere. 


Rectangular  Double  Cantilever  Specimen 

The  geometry  of  one-half  of  the  specimen  is  shown  in  Fig.  3  along  with 

the  details  of  the  finite  element  mesh  used  in  the  present  analysis.  Only 

2 

the  static  material  properties  E  =  3380  MN/m  and  v  =  0.33  for  the  present 
Araldite  B  material  are  used  in  the  present  analysis,  partly  due  to  the 
reason  of  comparing  the  present  results  with  those  of  Kobayashi  [9],  The 
presently  used  boundary  conditions  are  similar  to  those  in  [9].  The  experi¬ 
mental  data  for  crack-length  (£)  versus  time  history  that  is  simulated  in 
the  present  analysis  corresponds  to  specimen  No.  4  reported  by  Kalthoff 
et  al  [8],  In  reference  [8],  the  crack-initiation  stress-intensity  factor, 

Kj  ,  which  was  computed  fom  the  experimentally  measured  deflection  26  at  the 

3/2 

loading  points  by  using  the  formula  of  Kanninen  [17],  was  quoted  as  2.32  MN/m 
However,  the  value  of  6  was  not  quoted  in  [8].  By  using  the  formula  in  [17], 
the  load  point  deflection  to  create  an  initiation  K-factor  of  2.32  is  calculated 
to  be  5^.  However,  in  the  present  finite  element  method,  a  load  point  deflection 
of  (1.1)  was  found  necessary  to  impose  the  value  of  2.32  HN/m^  ^  on  the 
model.  This  discrepancy  may  be  due  to  the  approximate  nature  of  the  analysis 
in  [17]. 


The  crack-length  versus  time  history  (E  vs.  t)  for  specimen  ilo.  4  reported 
in  [8]  shown  is  Fig.  4.  The  crack-tip-velocity  versus  time  [E  vs.  t]  curve, 
obtained  by  numerically  dif ferentcating  the  (E  vs.  t)  curve  is  also  shown  in 
Fig.  4,  and  is  seen  to  be  in  good  agreement  with  that  shown  in  Fig.  6  of  [8]. 
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Kalthoff  et  al  [8]  measure  the  dynamic  stress-intensity  factors  by  applying 
the  method  of  caustics.  However  in  [8],  the  relation  between  the  diameter  D 
of  the  caustic  and  the  stress-intensity  factor,  is  based  on  the  static  soluticn 
for  the  stress-strain  field  around  the  tip  of  the  crack.  However,  Beinert 
et  al  [18]  show  that  if  the  above  relation  is  derived  from  a  dynamic  solution 
for  the  stress-strain  field  around  the  crack,  the  thus  computed  K-factors 
may  defer  from  those  in  [8]  by  a  factor  [1/F],  It  is  noted  that  only  the  basic 
singular  dynamic  plane  stress  solution  is  used  in  developing  the  relation 
between  and  D.  It  was  estimated  in  [18]  that  F  =  1.095  for  v  =  500  m/s,  and  F 
1.028  for  v  =  300  m/s,  where  v  is  crack-tip  velocity.  Further  it  is  noted  that 
the  epoxy  resin  Araldite  B  is  a  viscoelastic  material  whose  properties  are 
rate-dependent.  However,  in  the  present  analysis  as  well  as  in  [18]  the 
material  is  modeled  as  linear  elastic  through  out. 

It  is  seen  from  Fig.  4  that  the  presently  computed  dynamic  K-f actor,  Kd, 
agrees  well  with  the  experimental  data  [roughly  to  within  the  above  discussed 
factor  (1/F)]  until  the  time  when  the  crack-tip  begins  to  decelerate  as  per 
the  t( t)  curve  (which  is  obtained  by  numerical  differentiation  of  the  the 
experimental  data  for  E(t)).  In  the  time  interval  between  the  beginning  of 
crack-tip  deceleration  and  its'  eventual  arrest,  significant  differences  are 
noted  between  the  present  results  for  and  those  reported  in  [8]. 

To  understand  the  above  descrepancies  in  the  results  of  Fig.  4,  and  to 
study  the  effect  of  small  changes  in  the  assumed  Z  versus  t  curve,  the  analysis 
was  repeated  for  three  different  cases  of  assumed  £(t)  curves  as  shown  in 
Fig.  (5).  It  is  seen  that  the  curves  in  Fig.  (5)  represent  minor  deviations 

from  the  "experimental"  data  reported  in  f8].  It  is  also  noted  that  the  Z(t) 
curve  marked  as  "Data  3"  is  identical  to  the  one  used  by  Kobayashi  et  al  [h] 
in  their  "generation"  [9]  calculation.  The  presently  computed  curves, 


for  each  of  the  three  assumed  Data  for  £(t)  in  Fig.  (5),  are  shown  in  Fig.  (6). 
Also  shown  in  Fig.  (6)  are  the  £(t)  curves  corresponding  to  each  of  the  assumed 
E(t)  curves,  with  the  curve  for  'Data  3'  being  similar  to  that  in  [9].  The 
results  for  K^  computed  by  Kobayashi  in  his  "generation"  calculation  [9]  and 
the  experimental  results  of  [8]  are  also  shown  in  Fig.  (6).  It  is  seen  that 
significant  differences  exist  between  the  present  results  for  each  of  the 
three  E(t)  data  cases,  and  those  in  [8],  during  roughly  the  last  third  of  the 
crack-propagation  history.  Each  of  the  three  curves,  for  present  results, 
shown  in  Fig.  (6)  exhibit  a  pronounced  maximum  during  the  later  third  of 
crack-propagation  history.  The  results  of  Kobayashi  [9]  also  exhibit  a 
pronounced  maximum  for  K^,  which  however  is  seen  to  occur  earlier  than  in 
each  of  the  present  three  cases. 

Such  differences,  compared  to  the  experimental  data,  as  in  Fig.  6,  were 
also  noted  even  in  the  "propagation"  calculations  (in  the  sense  defined  in 
[9])  by  Kobayashi  [9],  and  Popelar  et  al  [10]. 

It  is  seen  from  Fig.  (6)  that  the  maximum  value  of  the  presently  computed 

Kj  is  the  largest  for  "Data  3"  while  it  is  smallest  for  Data  1.  Note  that  the 

rate  of  crack-tip  deceleration  is  the  most  severe  for  Data  3,  while  for  Data 

1,  the  deceleration  is  zero,  ie.,  the  crack-tip  is  still  propagating  with  a 

constant  velocity  of  295  m/sec.  Thus,  the  more  severe  the  crack-tip  deceleration 

is,  the  more  higher  is  the  maximum  in  the  computed  K^.  To  understand  the 

reasons  for  the  peak  in  K,  for  Data  1,  the  times  for  various  waves,  reflected 

d 

from  the  boundaries,  to  interact  with  the  propagating  crack-tip  are  computed 
and  shown  in  Fig.  (6).  With  A,  B,  and  C  denoting  the  three  boundaries  as 
marked  in  Fig.  3,  D^,  Dg  and  Dc  are  the  times  when  the  dilatational  waves 
reflected  from  the  boundaries  A,  B,  and  C,  respectively,  interact  with  the 
propagating  crack-tip.  Likewise,  and  Sg  are  the  times  for  the  shear  waves 


reflected  from  boundaries  A  and  B,  respectively,  to  interact  with  the  moving 
crack-tip.  It  is  interesting  to  note  that  the  computed  K^,  for  the  Data  1 
case,  begins  to  peak  at  the  instant  S^.  Further,  the  present  analysis  proce¬ 
dure  has  been  found  to  yield  excellent  correlations  in  several  constant  as  well 
as  non-constant  velocity  propagation  problems,  such  as  those  of  Broberg, 

Freund,  and  Nillson,  in  [6].  Thus  it  appears  reasonable  to  conclude  that 
the  results  in  Fig.  6  for  Data  1  may  be  accurate. 

To  further  understand  the  validity  of  the  results  for  the  most  severe  case 
of  Data  3,  the  variations  of  strain  energy  (U) ,  the  kinetic  energy  (T) ,  the 
fracture  energy  (F)  and  the  total  input  energy  (W)  were  studied  for  the  case 
of  Data  3  (which  as  seen  in  Fig.  6,  results  in  the  most  serious  descrepancy 


with  the  cited  experimental  results).  It  is  noted  that  in  the  present  analysis 
procedure,  the  dynamic  stress-intensity  factor  is  calculated  directly  as  a 
variable  in  the  finite  element  equations  [A, 5].  The  fracture  energy  is  computed 


Sj  =  [ 1- (v/Cj ) 2  ] 1/2 


In  Eq.  (8),  and  C2  are  the  dilatational  and  shear  wave  speeds  respectively. 

The  strain  energy  U  is  calculated  from  the  presently  computed  stress  and  strain 
data,  while  the  kinetic  energy  is  calculated  from  the  computed  data  for  velocities. 
Since  the  input  energy  to  the  specimen,  W,  is  a  constant  in  the  present  case, 
the  above  computed  F,  U,  and  T  should  add  up  to  a  constant.  It  is  seen  that 


the  error  in  (F  +  u  +  T)  as  compared  W  increases  almost  linearly  from  0.0/'  to 


about  4.5%  towards  the  end  of  the  computation.  Since  the  present  calculation 
was  carried  out  in  165  time-steps,  it  is  reasonable  to  presume  that  the  error 
in  each  time  step  is  thus,  roughly  0.0003%.  This  appears  to  give  enough 


credence  to  the  presently  computed  for  the  case  of  Data  3,  as  in  Fig.  6. 

The  results  in  Fig.  4  and  6  indicate  the  sensitivity  of  computed  K^,  in 
a  generation  calculation  (in  the  sense  defined  in  [9]),  to  the  input  data  for 
the  crack-tip  time  history,  I  vs.  t.  This  in  turn  points  to  the  extreme 
precision  with  which  the  history  E(t)  should  be  determined  in  an  experimental 
measurement.  Also,  the  validity  of  using  only  the  basic  singular  dynamic 
plane-stress  solution  in  developing  the  relation  between  K  and  the  caustic- 
diameter  as  in  [18]  appears  to  need  further  study. 

The  reason  for  the  descrepancies  in  the  presently  computed  results  and 
the  experimental  ones  for  Kd,  as  in  Figs.  (4)  and  (6)  may  be  explained,  in 
part,  by  the  fact  that  only  the  elastic  strain  energy,  kinetic  energy,  and 
fracture  energy  are  considered  in  an  energy  balance  relation,  of  the  type 
shown  in  Fig.  7,  in  the  present  procedure  based  on  a  linear-elastic  rate- 
independent  material.  However  for  viscoelastic  resins  of  the  type  of  Araldite 
B  used  in  the  experiments  [8],  the  viscous  dissipation  of  energy  may  also 
play  an  important  role  in  the  determination  of  K^.  This  may  distort  the 
comparison  of  such  experimental  results,  and  the  computational  results 
based  on  a  linear  elastic  material  behaviour. 

Part  II.  Thermally  Shocked  Cylinder  with  Meridional  Inner  Surface  Flaw 

Synopsis  of  the  Method  of  Approach:  Special  three-dimensional  hybrid 
crack  elements  are  used  to  model  the  immediate  vicinity  of  the  three- 
dimensional  crack  front,  and  the  conventional  20-noded  isoparametric  brick 
elements  are  used  to  model  the  remainder  of  the  cylinder.  The  special  crack 
elements  are  developed  through  a  hybrid  displacement  finite  element  procedure. 


i 


i 


Thus,  the  development  of  the  hybrid  crack  elements  is  based  on  a  three  field 

variational  principle  [11,12]  with  the  arbitrary  displacements  in  the  interior 

of  the  element,  a  displacement  field  at  the  boundary  of  the  element,  and  a 

Lagrange  multiplier  field  which  can  be  identified  as  the  traction  field  at 

the  boundary  of  the  element,  as  the  three  variables.  The  analytical  asymptotic 

solutions  for  the  displacement  field,  under  mixed  mode  (K^,  K^^.) 

conditions,  is  embedded  in  the  assumed  interior  displacement  field  of  the 

crack-element.  The  boundary  displacement  field  for  the  crack  element  is 

assumed  such  that  it  is  inherently  compatible  with  that  of  the  surrounding 

conventional  element.  The  third  field,  namely,  the  boundary  traction  field 

of  the  crack  element,  which  is  the  Lagrange  multiplier  to  enforce  the  equality 

of  the  boundary-value  of  the  assumed  interior  displacement  field  and  the 

independently  assumed  boundary-displacement  field  of  the  crack  element,  is 

assumed  such  that  it  includes  the  proper  (1/  r)  singularities  near  the  crack 

front.  From  the  details  of  the  present  finite  element  procedure  [11,12,13] 

it  can  be  seen  that  the  mixed  mode  stress- intensity  factors  (Kj,  K^j,  and 

at  various  points  along  the  crack  front  (which  are  denoted  by  a  master  vector 
* 

kg)  can  be  solved  for  directly  from  the  final  finite  element  equations: 

*  T  * 

~i  3  +  *2  bs  =  2l  (9) 

K2  q*  +  K3  k*  =  Q2  (10) 

where  (i=l ...  3)  are  the  corresponding  global  stiffnesses  and  (ot=l,2)  are 
the  corresponding  nodal  iorces. 

Results 

The  above  described  hybrid  crack  element  procedure  was  applied  to  analyse 

a  thermally  shocked  cylindrical  vessel,  of  commerical  geometry  with  outer  to 

inner  radii  ratio  of  (R  / R.)  =  (10/9),  and  containing  an  inner  surface  elliptical 

o  i 


flaw.  The  flaw  parameters  are:  (a/c)  =  0.2;  a/(RQ-R^)  =  0.6;  where  the 
parameters  are  defined  in  Fig.  8.  The  initial  temparature  of  the  vessel  is 
T^  and  inner  surface  of  the  vessel  is  assumed  to  be  instantaneously  reduced 
to  a  temparature  of  T^.  The  transient  temparature  distribution  in  the  cylinder 
is  obtained  from  [14]  as: 

t‘te  „ ;  -to2t  VWVWV™.)  , 

,e  "  ~T7, — ; — JZ — : — *wi7u  <u) 

i  E  m=l  J  (R.a  )  -  J  (R  a  )  o  i 

o  i  m  o  o  m 

where  U  (ra  )  =  J  (ra  )Y  (R  a  )  -  J  (R  a  )Y  (ra  ) ;  J  and  Y  are  Bessel  func- 
om  omoom  oom  om  o  o 

tions  of  the  first  and  second  kind,  respectively;  a  are  the  roots  of  U  (R.a  ) 

r  J  m  ..  o  i  m 

=  0;  K  is  the  thermal  diffusivity;  and  t  is  the  time.  For  purposes  of 

normalizing  the  final  solution  stress-intensity  factors,  the  maximum  value  of 

non-dimensional  hoop  stress  CS ^  =  [  (gA(|))mqvI/t'ta(T,-~TF)  I  is  taken  from  [14] 

to  be  5.1755  (where  E  and  a  are  the  modulus  of  elasticity,  and  the  coefficient 

of  thermal  expansion,  respectively). 

One  fourth  of  the  vessel  is  modeled,  and  the  finite  element  breakdown, 

which  is  shown  in  Fig.  9,  consists  of  380  elements  and  about  5600  degrees 

of  freedom.  The  solution  for  stress-intensity  factors  for  a  semielliptical 

surface  flaw,  with  (a/c)  =  0.2  and  a/(RQ-R^)  =  0.6,  is  shown  in  Fig.  10  along 

with  the  comparison  results  from  [14].  In  the  normalization  of  the  stress- 

intensity  factors,  the  normalizing  factor  0  ,  which  is  defined  earlier,  is 

used.  Thus  the  results  presented  in  Fig.  10  would  have  a  unit  of  stress. 

2 

Further,  the  ratio  (Kt/R^)  is  assumed  to  be  0.0001  in  the  present  problem. 

The  present  solution  is  seen  to  differ  from  that  in  [14]  by  about  10%  at 
the  free  surface  (0  =  0°)  of  the  flaw. 

However,  this  correlation  can  be  considered  to  be  good  for  purposes  of  practical 
application  of  an  engineering  fracture  theory.  From  the  present  results. 


as  well  as  the  additional  results  present  in  [13],  it  is  concluded  that  tin 


approximate  "alternating-technique"  solution  presented  in  [14]  may  be 
adequate  for  an  engineering  analysis  of  thermally  shocked,  flawed,  cylindrical 
pressure  vessels. 


Closure 

Novel  numerical  techniques  for  the  analysis  of  two-dimensional  fast 
fracture  situations,  and  three-dimensional  static  analysis  of  surface-flawed 
pressure  vessels,  have  been  presented  and  illustrated  through  solutions  of 
certain  important  problems  in  both  problem  areas.  The  present  numerical 
simulation  of  experimental  data  on  a  double-cantilever-beam  specimen  pointed 
to  the  need  for  highly  sophisticated  experimental  measurement  of  crack-tip 
velocity  history,  as  well  as  the  need  for  the  incorporation  of  more  realistic 
material  property  data  (for  experimentally  used  viscoelastic  resins  such  as 
Araldite  B)  in  a  computational  scheme.  Despite  the  tremondons  advances  made 
both  in  numerical  as  well  as  experimental  techniques,  much  needs  to  be  per¬ 
fected  in  both  the  areas! 
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Fig.  1  Finite  element  model  for  a  center-cracked  panel. 

Propagating  singularity  element  is  shown  my  hatched  markings. 

Fig.  2  Normalized  dynamic  stress  intensity  factor  variation  with  £(t). 

Fig.  3  Finite  Element  Model  for  a  Double-Cantilever-Beam  Specimen. 

Fig.  4  Normalized  dynamic  stress-intensity  factor  as  a  function  of  time, 
with  Kalthoffs  £(t)  curve  as  input. 

Fig.  5  Three  different  £(t)  curves  assumed  in  the  analysis. 

Fig.  6  Normalized  dynamic  stress-intensity  factor  as  a  function  of  time, 
for  each  of  the  three  assumed  £(t)  curves. 

Fig.  7  The  variation  with  time  of  different  energies  in  the  specimen. 

Fig.  8  Nomenclature  for  a  thermally  shocked,  internally  cracked  cylinder. 

Fig.  9  Finite  element  breakdown  for  the  crack  cylinder. 

Fig.  10  Stress- intensity  factor  variation  along  the  flaw-border  for  a 
thermally  shocked  cylinder. 
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